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ABSTRACT 

We use test-particle integrations to show that epicyclic motions excited by a pericentre 
passage of a dwarf galaxy could account for bulk vertical velocity streaming motions 
recently observed in the Galactic stellar disc near the Sun. We use fixed potential test- 
particle integrations to isolate the role of phase wrapping of epicyclic perturbations 
from bending and breathing waves or modes, which require self-gravity to oscillate. 
Perturbations from a fairly massive Sagittarius dwarf galaxy, Mg 2.5 X 10 10 M o , 
are required to account for the size of the observed streaming motions from its orbital 
pericentre approximately a Gyr ago. A previous passage of the dwarf through the 
Galactic disc approximately 2.2 Gyr ago (with a then more massive dwarf galaxy) 
is less effective. If phase wrapping of epicyclic perturbations is responsible for stellar 
streaming motions in the Galactic disc, then there should be variations in velocity 
gradients on scales of a few kpc in the vicinity of the Sun. 


Keywords: Galaxy: structure - Galaxy: kinematics 
and dynamics - Galaxy: disc. 


1 INTRODUCTION 

Significant bulk motions of stars, or streaming velocities, 
have recently been detected in large scale stellar surveys such 
as the RAdial Velocity Experiment (RAVE, Steinmetz et al. 
2006) and Large Area Multi-Object Spectroscopic Telescope 
(LAMOST, Cui et al. 2012) radial velocity surveys (Siebert 
et al. 2011; Williams et al. 2013; Carlin et al. 2013; Xu et 
al. 2015; Sun et al. 2015). Vertical wavelike structures in the 
stellar density distribution are also seen in the Sloan Digi¬ 
tal Sky Survey (SDSS) data (Widrow et al. 2012; Xu et al. 
2015). Curiously, the bulk motions of stars vary as a function 
of height above the plane, and the vertical bulk motions ex¬ 
hibit patterns of compression and rarefaction. The observed 
vertical velocity gradients could be due to heating from in¬ 
ternal perturbations, such as spiral arms or bars (Faure et 
al. 2014; Monari et al. 2015), or bending and breathing 
waves or modes (Widrow et al. 2012; Gomez et al. 2013) 
that could have been excited by a dwarf galaxy (Widrow 
et al. 2014). Vertical density and velocity structures arise 
in simulated Milky Way discs perturbed by a Sagittarius 
sized dwarf galaxy (Gomez et al. 2013) that could also have 


induced warped and ringed structures (such as the Mono- 
cerous Ring; Newberg et al. 2002) in the outer Galaxy (e.g., 
Kazantzidis et al. 2008; Younger et al. 2008; Quillen et al. 
2009; Purcell et al. 2011). 

Pericentre approaches and passages through the Galac¬ 
tic disc by the Sagittarius dwarf galaxy would have excited 
both vertical and radial epicyclic motions in stars. A per¬ 
turbation in the disc could affect a localised region where 
stars are pushed to the same phase in their epicycles. Over 
time, these perturbed stars would see a large spread in 
phases, called phase wrapping, due to the dependence of 
vertical and radial (epicyclic) oscillation frequencies on or¬ 
bital angular momentum, eccentricity and inclination (e.g., 
Minchev et al. 2009), giving spiral and warped structures in 
the disc (Quillen et al. 2009). However, the perturbation of a 
dwarf galaxy could also excite vertical bending and breath¬ 
ing waves or modes in the disc (Weinberg 1991; Gomez et al. 
2013; Widrow et al. 2014; Widrow V Bonner 2015). N-body 
simulations would be expected to give a more realistic sim¬ 
ulation of perturbations to the Galactic disc caused by the 
Sagittarius dwarf galaxy. However, with N-body simulations 
it is difficult to differentiate between excitation of bending or 
breathing waves or modes from phase wrapping of epicyclic 
motions. Waves and vibration modes are a property of a self- 
gravitating disc (Widrow et al. 2012), and so would not be 
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present in a non-interacting test-particle simulation where 
particle motions are integrated in a fixed Galactic potential. 
Here we use test particle integrations to isolate and study 
the role of epicyclic phase wrapping and explore which struc¬ 
tures in the local velocity field might be a result of vertical 
and radial epicyclic motions excited in the disc by a previous 
pericentre or a previous passage through the Galactic disc 
of the Sagittarius dwarf galaxy. 


2 TEST-PARTICLE INTEGRATIONS 

We integrate particle orbits for two different sets of initial 
conditions using a fixed gravitational potential for the Milky 
Way. Particle orbits are integrated using the python library 
galpy (Bovy 2015) using the MWPotential2014 model grav¬ 
itational potential. This potential is designed to give a re¬ 
alistic model for the Milky Way and is described in detail 
by Bovy (2015) in his section 3.5. Our procedure is as fol¬ 
lows: We generate a thin disc of 10 7 test particles. We per¬ 
turb the velocities of each particle, instantaneously approx¬ 
imating perturbations caused by an encounter from a dwarf 
galaxy. We then integrate the particles to the present day in 
the static Galaxy potential alone. Each set of integrations 
required 3-5 days of computations on a Mac-Mini computer 
with a 2.4 GHz IntelCore 2 Duo processor. The perturba¬ 
tions are applied instantaneously, rather than directly in¬ 
tegrated as a function of time, so as to separate between 
the excitations caused by the perturbations and subsequent 
phase wrapping in the fixed background Galactic potential. 

We consider the role of two perturbations from the 
dwarf galaxy, separately, so that we can contrast and com¬ 
pare their different roles in exciting epicyclic motions in the 
stars. We run two separate integrations, one beginning with 
perturbations caused by the passage of the dwarf galaxy 
through the Galactic plane approximately 2.2 Gyr ago, and 
the second integration, beginning with perturbations caused 
by the dwarf galaxy’s pericentre approximately 1.1 Gyr ago. 

We work in the right-handed Galactocentric coordinate 
system used by Law & Majewski (2010). The relation be¬ 
tween Galactic latitude and longitude (on the sky and with 
origin at the Sun) and the Cartesian Galactocentric coordi¬ 
nate system with origin at the Galactic Center is illustrated 
in Figure 1. This coordinate system gives a current location 
of the Sagittarius dwarf at X, Y, Z = (19.0, 2.7, —6.9) in kpc 
with a distance from the Sun of ~ 28 kpc and Galacto¬ 
centric solar radius Rq = 8.0 kpc. In this coordinate system 
the current location of the Sun is (— Rq , 0, 0) and this differs 
from that used by Carlin et al. (2013) and Sun et al. (2015) 
who place the Sun at (R©,0,0). We also use Galactocen¬ 
tric cylindrical coordinates (R, 0, Z), with Galactic rotation 
0 > 0, and in these the location of the Sun is (R©,0©,O) 
with 00 = 7r, differing from Carlin et al. (2013) and Sun et 
al. (2015) who adopt 0© = 0. We work in units of kpc, M© 
and velocities are given in km s -1 . The adopted circular ve¬ 
locity at Rq is Vq = 220 km s _1 and the period of Galactic 
rotation at Rq is P — Rq/Vq = 0.23 Gyr. 


North Galactic 



Figure 1 . The top panel shows Galactic latitude and longitude 
(on the sky). The bottom panel shows the Galactocentric carte¬ 
sian coordinate system. 


2.1 Thin disc prior to perturbation by the dwarf 
galaxy 

Prior to perturbation by a dwarf galaxy encounter, parti¬ 
cles are initially evenly distributed in a thin disc. Work¬ 
ing in Galactocentric cylindrical coordinates we start by 
choosing azimuthal angle, 0, from a uniform distribution 
0 ^ 0 ^ 27 t. Guiding radii are chosen to be uniformly dis¬ 
tributed 5 ^ R g ^ 30 kpc, giving a surface density pro¬ 
portional to 1/R. This has a shallower radial gradient than 
an exponential distribution but allows us to more efficiently 
study the structure of the outer Galaxy by increasing the 
number of test particles in that region. 

After the guiding radius R g is chosen, the initial radius, 
radial and tangential velocity components are chosen using 
a first order (in radial action variable) epicyclic approxima¬ 
tion. The epicyclic amplitudes a r are chosen to be uniformly 
distributed between zero and a rma£C . Epicyclic angles, 0 r , are 
chosen from a uniform distribution with 0 ^ (j) r ^ 2tv. Initial 
radii and radial velocities are set to 

R = a r sin (cf) r ) + R g 
Vr = a r K(R g ) cos(0 r ). (1) 

The vertical component of the orbit’s angular momentum is 
Lz = R g V c (R g ) and the initial tangential velocity compo¬ 
nent computed as Vq — Lz/R- Here V c {R g ) is the rotation 
curve and ft(R g ) is the radial epicyclic frequency, and these 
are computed using the galactic potential MWPotential2014 
using the vcirc and epifreq routines from galpy. The ini¬ 
tial vertical positions and velocities are similarly chosen by 
choosing an epicyclic amplitude, a z , uniformly distributed 
between zero and a zrn ax , and vertical epicyclic angle (f) z , uni¬ 
formly distributed between zero and 27r, then setting initial 
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Z and vertical velocity component Vz to be 
Z = a z sin(4> z ) 

Vz = a z v(R g ) cos(^). (2) 

Here v(R g ) is the vertical oscillation frequency, which is 
computed in the MWPotential2014 galactic potential using 
the verticalfreq routine from galpy. Setting a rmacc = 1.6 
kpc and a zrn ax =0.4 kpc, the resulting disc has velocity 
dispersions <jr ~ 20 km s _1 , <jz ~ 10 km s _1 , and a$ ~ 15 
km s _1 at R — Rq that are consistent with thin disc values 
computed in the Solar neighbourhood (comparing to Fig¬ 
ure 31 by Nordstrom et al. 2004). The parameters a r max 
and CLzmax are independent of galactocentric guiding radius, 
(R g ), corresponding to a thin disc with scale height indepen¬ 
dent of radius. As v/n varies from 2.0 at Rq to 1.56 at 25 
kpc the ratio az /ctr is 0.5 at Rq and drops to 0.4 at 25 kpc. 
The expected vertical density profile p(z) is logarithmically 
sensitive to z/a zrn ax with a cusp at the midplane. From the 
initial stellar distribution, we numerically measured the ver¬ 
tical dispersion of the initial density profile at Rq finding a 
width yj (z 2 ) = 163 pc. We also numerically measured the 
asymmetric drift v a = (V c — Vq) ~ 3 km s _1 at Rq. 

2.2 Perturbations from a Dwarf Galaxy 

The orbit of the Sagittarius dwarf galaxy has been inferred 
from modeling the continuous stream of debris in the Sagit¬ 
tarius dwarf galaxy’s tidal tails (e.g., Johnston et al. 2005; 
Law X Majewski 2010; Purcell et al. 2011) and resembles a 
trefoil knot (see Figure 1 by Johnston et al. 2005 and our 
illustration in Figure 2). The Sagittarius dwarf galaxy nu¬ 
cleus is currently near pericentre in its orbit. The Sagittarius 
dwarf galaxy nucleus previously passed through the Galactic 
plane about 0.4 Gyr ago, but at large Galactocentric radius 
R > 50 kpc. The previous pericentre occurred about 1 Gyr 
ago and at positive Z in our Galactocentric coordinate sys¬ 
tem (see lower left panel of Figure 8 by Law & Majewski 
2010). Our illustration in Figure 2 labels the last pericentre 
as El. About 2 Gyr ago, the dwarf galaxy’s nucleus passed 
through the Galactic plane when it was nearly at pericentre. 
This event is labelled as E2 in our illustration. In Figure 2 
we plot a model orbit by Chakrabarti et al. (2014) (that la¬ 
belled E in their Table 1), that matches the observed proper 
motions for stars in the Sagittarius dwarf galaxy nucleus. Us¬ 
ing this orbit, we list in Table 1 the positions and velocities 
of the dwarf galaxy nucleus at the El pericentre and the E2 
passage through the Galactic plane. These two events could 
have caused large perturbations to the disc stars and we use 
them to generate velocity perturbations for our test-particle 
integrations. 

After generating a thin disc we use an instantaneous 
hyperbolic orbit approximation (often used to derive dy¬ 
namical friction and gravitational heating rates, see Binney 
X Tremaine 1987 section 7.1) to perturb the velocities of 
each particle. Two sets of integrations are done, each using 
a separate encounter listed in Table 1. We use a hyperbolic 
orbit to compute the velocity perturbations on each parti¬ 
cle (applied instantaneously at the beginning of the orbital 
integration). Our procedure for computing the velocity per¬ 
turbation from the dwarf encounter positions and velocities 
in Table 1 is described in more detail in our appendix A. 



Figure 2. Illustrating the Sagittarius dwarf orbit from the E 
model by Chakrabarti et al. (2014). The current position of the 
dwarf in its orbit is labelled ” SagD” and the passage through the 
galactic plane 2.2 Gyr ago labelled with E2 and a black arrow. 
The previous pericentre, 1.1 Gyr ago, is labelled El and shown 
with another black arrow. Axes are in Galactocentric coordinates 
and in units of kpc. 

Using this approximation the velocity perturbation Av is 
proportional to the mass of the dwarf galaxy, M^. Hence, 
perturbations from larger or smaller dwarf galaxy masses 
can be estimated by scaling the resulting velocity perturba¬ 
tions. 

The progenitor dwarf galaxy mass for the E model by 
Chakrabarti et al. (2014) is Md = 10 lo M© and lighter than 
both the ‘light’ (3 X 1O 1O M 0 ) and ‘heavy’ (10 n M©) Sagit¬ 
tarius dwarf galaxy masses explored by Purcell et al. (2011). 
Purcell et al. (2011); Chakrabarti et al. (2014) take into ac¬ 
count dynamical friction and tidal stripping of the dwarf 
galaxy. The dwarf galaxy mass is similar to its initial mass 
until pericentre near E2, approximately 2 Gyr ago, at which 
time the dwarf loses about half of its mass (see Figure S2 
by Purcell et al. 2011). The dwarf galaxy mass then remains 
constant until the El pericentre. Consequently, we take the 
mass of the dwarf to be equal to its progenitor mass for the 
E2 encounter and half that at the El encounter. For the E2 
encounter we use a dwarf galaxy mass of Md E2 = 5 x 10 10 M© 
(which implies Md E1 = 2.5 x 10 lo M©), in between the light 
and heavy models by Purcell et al. (2011), but heavier than 
the E model by Chakrabarti et al. (2014). 

2.3 Initial mean velocities 

We generate 10 million particles for the El and the E2 en¬ 
counters, separately. Following generation of the thin disc 
and the instantaneous velocity perturbation due to the El or 
E2 encounters, (and prior to orbit integration) we compute 
mean velocities in 0.25 kpc square bins in X and Y, summing 
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Table 1. Parameters for the two dwarf galaxy encounters labelled in Figure 2 


Encounter 

t (Gyr ago) 

Mass (xlO 10 Mq) 

v (km s 1 ) 

r (kpc) 

El 

1.1 

2.5 

(-339, -44, 76) 

(4,8,18) 

E2 

2.2 

5.0 

(4,-111,-337) 

(-21,-4,0) 


together the velocity components for all particles within each 
bin. Figure 3 shows the initial mean radial velocity compo¬ 
nent, (Vr), the mean vertical velocity component (Vz}, and 
the mean tangential velocity component subtracted by the 
circular velocity as a function of radius, (Vo) — V C irc(R ), for 
the El and E2 encounters. The E2 encounter has velocity 
vector nearly but not exactly perpendicular to the Galaxy 
midplane. The velocity vector is tilted slightly so that the 
orbit is at positive Z above the midplane at negative Y and 
vice-versa at positive Y. The tilt of the dwarf’s orbit gives a 
negative Vz perturbation to the particles in the disc at pos¬ 
itive Y and vice-versa at negative Y. We see this as nearby 
red and blue regions in the (Vz) panel in Figure 3. These 
nearby particles differ in vertical epicycle angle, <p z , by 7r 
(see equations 2). The situation is different for the El en¬ 
counter. Because the dwarf galaxy passes above the Galactic 
plane in this encounter, both (Vr) and (Vz) are in the same 
direction, though (Vo) has both positive and negative pertur¬ 
bations (see top panels in Figure 3). Because (Vz) is positive 
the perturbed particles have the same phase <\> z ~ 0 and the 
size of (Vz) gives vertical epicyclic amplitude a z ~ 0.05 in 
units of radius (using a typical Vz and equations 2). The sign 
change in the Vo perturbation gives epicyclic angle (f> r ~ tt/2 
where (Vo) is lower than the circular velocity and (f> r ~ — tt/2 
where (Vo) exceeds the circular velocity (using equation 1) 
and epicyclic amplitude a r ~ 0.15 in units of radius. Even 
though the dwarf mass used for the El encounter is half of 
that in the E2 encounter, the Vz perturbation is larger. This 
is because the El encounter has its closest approach above 
rather than in the Galactic plane. Even though the dwarf 
mass is lower for the El encounter, it is more effective at 
exciting vertical epicyclic motions. 


3 VELOCITY DISTRIBUTIONS AFTER 
ORBITAL INTEGRATION 

Using galpy’s Dormand-Prince integrator, the initial parti¬ 
cle distributions (with mean velocities shown in Figure 3) 
are then integrated to the present time (2.2 Gyr for the E2 
integration and 1.1 Gyr for the El integration). In Figures 

4 and 5 we show mean velocities, similar to those presented 
in Figure 3, as well as the particle density distributions as a 
function of A, V for both El, and E2 encounters, but after 
orbit integration to the present time. The density distribu¬ 
tions shown in Figures 4 and 5 resemble those illustrated 
by previous works using particle integrations (Quillen et al. 
2009) and N-body simulations (Purcell et al. 2011; Gomez et 
al. 2013; Widrow et al. 2014). We confirm that trailing spiral 
structure or overdensities can be induced by dwarf galaxy 
encounters, if the dwarf galaxy mass is similar to 10 10 Mq. 
Because self-gravity is not present in our simulations, the 
spiral structure is due to the radial epicyclic motions excited 


by the encounters that are evident in the initial perturba¬ 
tions to Vr and Vo shown in Figure 3. 

In Figures 4 and 5 we also show the mean Z value as a 
function of X, Y. This is computed in 0.25 kpc square bins 
in the same way that we compute the mean velocity com¬ 
ponents. The (Z) subpanels in Figures 4 and 5 show that 
there are large regions of the disc with mean height above or 
below the galactic midplane, implying that the galaxy disc 
has become warped. This phenomenon has also been seen 
in previous simulations (e.g., see Figure 5 by Gomez et al. 
2013). The warp arises as vertical amplitudes induced by the 
encounter (and seen in the (Vz) distributions just after the 
encounters in Figure 3) have been sheared due to differential 
precession and the radial gradient of the vertical epicyclic 
frequency v. We note that this warp is not self-consistent 
as we use a static galactic potential, whereas a warp in the 
disc would influence the potential. The El encounter seen 
in Figure 4 shows ring - like structure in (Z) that is found 
in similar regions as and resembles the Monoceros and Tri¬ 
angulum - Andromeda overdensities as observed by Xu et 
al. (2015) and Price-Whelan et al. (2015). We note that the 
warp for E2 encounter appears to be more twisted than that 
of the El integration. The E2 encounter occurs 2.2 Gyr ago 
and there are more rotation periods for the phase wrapping 
to take place. 

3.1 Breathing and Bending Coefficients 

Breathing and bending modes arise from the self-gravity of 
the Galactic disc. Breathing modes signify that stars above 
and below the midplane move in opposing vertical direc¬ 
tions, while bending modes imply stars above and below the 
midplane move en masse upwards or downwards (Sun et al. 
2015). To quantify bending and breathing modes Widrow et 
al. (2014) fit a function to the mean vertical velocity as a 
function of Z in each A, Y bin, 

(Vz){x,y,z) = A z (x,y)z + B z (x,y) (3) 

(their equation 24). Here Bz(x,y) is not equivalent to our 
computed (Vz) as it is equal to the mean value of Vz at 
Z = 0, and the density distribution may not be symmetric 
about Z — 0. We use a linear regression to fit for the Az 
coefficient in each 0.25 kpc square bin in A, Y. Widrow et al. 
(2014) use Az to describe breathing modes and Bz to de¬ 
scribe bending modes. The Az and Bz coefficients are also 
shown in subpanels of Figures 4 and 5. We have measured 
Az and Bz coefficients with sizes of several km s -1 kpc and 
km s _1 respectively for both encounters. Figure 4 shows 
that the structure of the Az and Bz coefficients in our El 
encounter resembles the spiral structure seen in Figures 10 
and 12 in Widrow et al. (2014), though it is necessary to 
note that these figures from Widrow et al. (2014) occur at 
different times from our encounters after the initial pertur¬ 
bation. As we have carried out non-interacting test-particle 
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Figure 3. Mean particle velocity distributions as a function of X , Y in the Galaxy prior to orbital integration. The initial thin discs have 
been instantaneously perturbed using an instantaneous hyperbolic orbit approximation (described in appendix A). In the top panels, 
we show the mean velocities induced by the pericentre encounter 1 Gyr ago (El encounter). The top left panel shows the mean radial 
velocity component, the top middle panel, the mean vertical velocity component and the top right panel the mean tangential velocity 
component subtracted by V c i rc {R). Each panel displays the position and direction of the Sagittarius Dwarf in X , Y with a black dot and 
an arrow. For scale the current position of the Sun is also shown, however the Sun would have been at a different azimuthal angle during 
the encounter. Galactic rotation is counter-clockwise. On the bottom a similar set of panels show mean velocities but for the passage 
through the Galactic plane 2 Gyr ago (E2 encounter). Parameters for the two encounters are listed in Table 1. Nearby particles differ in 
vertical epicyclic angles by 7r for the E2 encounter, whereas they all have the same phase in the El encounter. Even though the dwarf 
mass is lower for the El encounter, it is more effective at exciting vertical motions. 


integrations in a fixed potential, these coefficient values and 
structure cannot be due to bending or breathing waves or 
modes. Rather, they can only be due to phase wrapping of 
epicyclic amplitude perturbations. 


The perturbing dwarf galaxy is more massive for the E2 
encounter, however there is longer time for phase wrapping 
to occur and the perturbed motions vary in phase. Values for 
the Az coefficient in Figure 5 range from -5 to 5 km s -1 kpc 
but much of the substructure is erased because the structure 
is so tightly wound. For the El simulation, structure in Az 
is more coherent and ranges from -10 to 10 km s -1 kpc. The 
extent of variations seen in Az in our integrations is nearly 
the same size as those observed locally; Carlin et al. (2013) 
and Sun et al. (2015) measure Az values in the range -10 
to 10 km s -1 kpc (see Figure 13 by Sun et al. 2015 and 
Figure 3 by Carlin et al. 2013). Variations in the value of 
(Vz) are approximately ±10 km s -1 for both simulations 
and similar to those observed (see Figure 14 by Williams et 
al. 2013). While phase wrapping from a relatively massive 
dwarf galaxy with a single encounter can cause measurable 
values for the Az coefficient, it may not account for the full 
extent of the observed values. 


3.2 The vertical gradient of the mean radial 
velocity component 

To quantify the sensitivity of the radial velocities with height 
above the Galactic plane we also fit a linear function to the 
mean radial velocities 

C Vn)(x,y,z) = A R (x,y)z + B R (x,y) 

similar to equation 3 but using Vr and with coefficients 
Ar , Br. These maps are computed in the same way as for 
the Az,Bz coefficients and also shown in Figures 4 and 5. 
We can see from these Figures that there are strong corre¬ 
lations between the maps. 


3.3 Dynamical Interpretation 

Our original velocity perturbations (shown in the initial ve¬ 
locity perturbations in Figure 3) are localised near a par¬ 
ticular azimuthal angle (near closest approach) that we can 
denote O p , however, Vr and Vz velocity perturbations also 
extend over a range of radius. We can roughly group stars 
into two sets: those that are perturbed, and those that re¬ 
main in circular planar orbits. After 1 or 2 Gyr (depending 
on the encounter) the dependence of the angular rotation 
rate on radius, Q(R) shears the distribution of perturbed 
stars. Meanwhile the stars oscillate vertically and radially, 
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X (kpc) 


X (kpc) 


X (kpc) 


Figure 4. Distribution of mean velocities and Az , Bz , Ar, Br parameters as a function of X,Y for the El encounter after orbit 
integration to the present time. Top row from left to right, Az,Bz and (Vz), middle row from left to right Ar,Br and (Vr), bottom 
row, (Z), (Vq) — V c i rc (R) and difference between current and initial number densities. Velocities and Bz,Br are in km s -1 , and the 
Az , Ar coefficients are in km s —1 kpc. The difference between current and initial densities is shown in numbers of particles per 0.25 kpc 
square bin. 


at frequencies that depend on their mean radius or angular 
momentum. 

For an initial perturbation at 0 P , after a time AT the 
perturbation will be wrapped in 27r in azimuthal angle across 
a radial distance 


A R ~ 


27T 

Q, r A T 



(4) 


where 0, r is the radial derivative of the angular rotation 
rate and P ~ 0.23 Gyr is the rotation period at the Sun’s 
Galactocentric radius and we have approximated the rota¬ 
tion curve as a flat one with V C (R) = constant. For our per¬ 
turbation with AT = 1 Gyr ago, A R ~ 1.8 kpc. That means 
that every 1.8 kpc in radius we should encounter stars that 
were initially in the peak of the perturbed region, however 
this neglects radial motions in the perturbed stars. 

The initial perturbation region, in cylindrical coordi¬ 
nates, is illustrated in Figure 6 in the top left panel as a 


grey bar. Ignoring the radial motions of the perturbed stars, 
after a time AT later, the distribution of perturbed stars is 
tilted, as shown in the top panel second from left. Stars out¬ 
side the perturbed region have little radial motion, so they 
remain in the region outside the tilted bar but they shear 
at the same speed. For the El encounter, initial Vz > 0, 
giving the perturbed stars an initial epicyclic phase 4> z « 0. 
The grey bar shown in Figure 6 will oscillate vertically as 
it tilts, after a time AT the perturbed stars are at differ¬ 
ent heights (shown in top panel second from right). Peak to 
peak vertical maxima should occur across a radial distance 


ar " = r ItV 


(5) 


and computing this distance for AT = 1 Gyr and v/VL ~ 2.6 
(using galpy’s Milky Way potential MWPotential2014 ) we 
estimate A R v ~ 0.75 kpc. This we expect is the radial dis¬ 
tance between two maxima in Z in the inner region of the 
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Figure 5. Similar to Figure 4 except for the E2 encounter. 
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galaxy, and it is approximately consistent with what is seen 
near Rq in the lower left panel of Figure 4 showing (Z). 
We compare the peak to peak radial distance from epicyclic 
phase wrapping to the wavelength predicted for bending 
waves, A R v ~ 10 kpc (Weinberg 1991, see Section 6). The 
radial distance between peaks in Z can be much smaller than 
the wavelength predicted for bending waves. 

The perturbed stars also oscillate radially. If there is a 
single phase for (f> r for the perturbed stars, then the distri¬ 
bution of the perturbed stars wiggles, as shown by the green 
curve on the top right panel. The wiggles in the green curve 
can lie in the same regions as stars on circular orbits (yellow 
background). There are overlap regions that would contain 
higher surface densities of stars accounting for the similar 
morphology of the density panel and velocity panels in Fig¬ 
ure 4. For some R , 0 there are two populations of stars: stars 
in circular orbits and stars that were perturbed by the en¬ 
counter. As the stars perturbed by the encounter can reach 
larger distances above and below the Galactic plane, the 
combined population can exhibit a vertical velocity gradi¬ 
ent. Likewise, there are regions where perturbed stars out¬ 


number unperturbed stars, leading to underdensities that 
can still exhibit velocity gradients. 


The El encounter also gives velocity perturbations in 
Vr and Vq , resulting in a region of the galaxy with stars that 
all have an epicyclic angle (f> r ~ —7t/ 2 for 0 slightly lower 
than the center 0 P and (j) r ~ tt/2 for 0 slightly higher than 
the center, the change in sign resulting from the azimuthal 
variation in Vq (see Figure 3). We can consider a localised (in 
0) increase in epicyclic amplitudes a r and a z , with (j) z m 0, 
near the center and <j) r — — 7t/ 2 but flipping sign at 0 P . This 
situation is represented in the lower left panel in Figure 6 
and after shearing in 0 in the lower panel, second from left. 
Vertical oscillations are in phase and shown in the lower 
panel second from right. Because the epicyclic frequency, 
(k) and vertical oscillation frequency (i/) depend on radius 
as the region of perturbation is stretched, it also oscillates 
radially and up and down. A perturbation in the radial di¬ 
rection and initially all in phase in 0 r will be wrapped in its 
radial epicyclic by 2tt (corresponding to two particles each 
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Figure 6. A perturbation localised in © but extending in R is shown on the top left panel. After a time AT the perturbation shears 
due to differential rotation (top panel second from left). Taking into account the frequency of vertical oscillations, the perturbed stars 
move up and down (top panel second from right). Taking into account a perturbation in radial epicyclic and assuming the perturbation 
has a single phase in radial epicyclic angle </> r , the radial motions deform the distribution of perturbed stars giving the green curve on 
the top right panel. Where the green curve lies outside the white bar, perturbed stars are found at the same R , © as unperturbed stars 
and the combined populations can give vertical gradients in the velocity components. The El encounter initially has both positive and 
negative changes to Vq giving epicyclic phases that are both positive and negative (shown on the bottom left panel) and after shearing 
due to differential rotation (shown in the bottom panel second from left). Vertical oscillations are in phase and shown in bottom panel 
second from right. The radial excursions from each bar have opposite sign and so can overlap one another as well as unperturbed stars 
(shown on bottom right panel). 


at pericentre) across a radial distance 


ar ~= r ^t-. 


( 6 ) 


Using k/Q ~ 1.4 we estimate A R K ~ 1.3 kpc near the Sun. If 
all perturbed particles were initially in phase then we would 
expect the distance between radial oscillation maxima to 
be 1.3 kpc, however the <fi r varies by almost i r in the initial 
perturbation and this gives about twice as many oscillations. 
Furthermore, perturbed stars can overlay each other as well 
as overlay unperturbed stars (as shown in bottom right panel 
in Figure 6). The two effects likely account for the increased 
number of red and blue regions in the Ar and Br panels 
compared to the Az and Bz panels shown in Figure 4. 

In Figure 4 the (Vz) and (Z) subpanels resemble each 
other except they are 90° out of phase, as would be expected 
from vertical oscillations of the perturbed stars. The Az 
slope coefficient measures the gradient of (Vz) with Z. In 
regions where stars in the midplane coexist with (or overlap) 
perturbed stars, a gradient is measured when the perturbed 
stars lie above or below the Galactic plane and there are 
unperturbed stars in the midplane. Where the perturbed 
particles are above the plane (Z > 0) but moving downward 
(Vz < 0), we estimate the gradient from these stars and 
the unperturbed ones in the midplane with Vz ^ 0, giving 
approximately Az ~ ( Vz) I(Z) and Az is negative. For stars 
below the plane (Z < 0) with Vz < 0, Az > 0. Thus, we 
expect the Az map contains twice as many color features as 
the (Z) and (Vz) maps and indeed that is what we see on 
the top left in the Az panel of Figure 4, when compared to 
the (Vz) and (Z) panels. 

Similar phenomena are seen in Figure 5, showing the 
E2 encounter, except there are many more features in each 


panel. This is expected as the time to wrap in phase is twice 
as long, but in addition the initial perturbation contains re¬ 
gions differing in phase (j) z and (j) r also contains additional 
structure as initial Vr and Vz distributions are more com¬ 
plex than for the El encounter (see Figure 3). Because the 
E2 encounter is less effective at inducing vertical oscilla¬ 
tions and much of induced structure is erased due to the 
long timescale and higher complexity of the epicyclic phases, 
hereafter we primarily discuss the El encounter. 

Noticeable correlations between the density of stars and 
different velocity components exist in our simulations. Fig¬ 
ures 4 and 5 show that regions that have extreme velocity 
means are the same as regions with extreme gradients, and 
these are the same as those with high or low surface den¬ 
sities. Perturbed stars lie chiefly in regions where the mean 
velocity component is nonzero and there are overdensities 
or underdensities. As perturbed stars lack circular orbits, 
they have radial velocity components that affect their tan¬ 
gential velocities; in general, positive (Vr) gives negative 
(Vo) — Vdrc(R), and this occurs where there are overdensi¬ 
ties in our integrations. Differences in bulk radial velocities 
can result in separated overdensity streams, where bands 
of stars share the same density but have split off from each 
other, such as those seen around ( X , Y) — (—15,15) kpc and 
(. X , Y) s=s (5,10) kpc in the density panel in Figure 4. These 
disparities in (Vr) also produce bifurcations, or ‘forked - 
tongue’ patterns, in (Vz), seen in the upper right panel in 
Figure 4. 

Negative values of (Vo) — V C i rc (R) (corresponding to pos¬ 
itive asymmetric drift) are expected from a steady state disk, 
though here the density, mean velocity components and dis¬ 
persions are time dependent. Because the initial density dis- 
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tribution is proportional to 1/R we expect that eccentric 
particles are more likely to affect the velocity mean near 
apocenter where Vq is lower than that of particles in cir¬ 
cular orbits, rather than near pericenter where Vq is higher 
than that of particles in circular orbits. 

Vertical velocity gradients are seen in our numerical in¬ 
tegrations and they can be explained by vertical and radial 
epicycle motions excited by the encounters. The vertical gra¬ 
dients arise in regions where perturbed and unperturbed 
stars (or different source regions of perturbed populations) 
lie in the same X , Y region. Because they arise due to overlap 
of different populations we expect (and see in our integra¬ 
tions) strong correlations between the distributions of the 
different velocity components and the density. The number 
of peaks in (Z) or (Vz) can roughly be estimated using a 
shearing timescale based on the oscillation frequency and 
the timescale since the perturbing encounter. 


3.4 Local Velocity Gradients 

To mimic observations of nearby stars in the Galaxy we can 
extract test particles in small regions of the Galaxy. Sun et 
al. (2015) (in their Figures 13 and 14) showed linear fits to 
{Vz) with Z at different locations in X , Y near the Sun. We 
construct here a similar figure to illustrate that our integra¬ 
tions display gradients in the vertical velocity component. 
A matrix containing numbers of particles is constructed in 
(8 ^ R ^ 10 kpc, 176° ^ 0 ^ 184°, — 1 ^ Z ^ 1 kpc) with 
binning size (AR, A0, A Z) = (0.25 kpc, 2°, 0.25 kpc). After 
generating the number density matrix, we apply the same 
outlier rejection criteria used by Sun et al. (2015), namely 
requiring —200 < Vr < 200 km s -1 , 0 < Vq < 400 km s -1 
and —200 < Vz < 200 km s -1 to compute a velocity 
matrix in (Vz)] we then divide the velocity matrix by the 
number density matrix to yield a bulk or mean Vz matrix 
in R , 0 and Z. We use linear regression to find the best fit 
line of 

(Vz) (R, 0, z ) = A Z (R, 0)z + B Z {R , 0) 

for each pair of (R, 0). The individual values for {Vz){z) 
are shown in Figure 7 at © = 180° along with the linear fits 
giving us Az and Bz coefficients; the A, Y axes of each sub¬ 
panel correspond to — 2 ^ Z ^ 2 kpc and — 30 ^ (Vz) ^ 30 
km s _ 1 , respectively. In the previous section we attributed 
the velocity gradients to overlap of populations of perturbed 
and unperturbed stars. However, most of the panels in Fig¬ 
ure 7 show smooth variations in mean velocities as a function 
of height above and below the Galactic midplane. Had we 
started with a planar thin disc in circular orbits we would 
have seen folds in the disc and these would not have given 
smooth velocity gradients. We can attribute the smooth 
slopes to the velocity dispersion in our initial disc. Sun et 
al. (2015) (and others) measured different velocity gradients 
in different regions of the galaxy. Figure 7 shows that there 
are rapid variations in velocity gradients with position in 
our simulated Galaxy. Hence some or all of the fine struc¬ 
ture measured by Sun et al. (2015) and Carlin et al. (2013) 
might be real, and we might expect that future observa¬ 
tions extending the number of stars, precision and distance 
of stars will uncover even more structure in the velocity held. 

We now explore whether our integration can exhibit 


the patterns seen as a function of Z in recent observational 
studies. A key feature seen in the local velocity distribu¬ 
tions is a strong vertical gradient in (Vz) with negative Vz 
at positive Z and a slope of Az ~ —5 to —10 km s -1 kpc 
(Carlin et al. 2013; Williams et al. 2013). Regretfully, in 
Figure 4 Az > 0 at (—R©, 0,0) and so is not consistent 
with the observed pattern of gradients near the Sun. As the 
gradients are comprised of tightly wound structures, they 
are sensitive to the time since the encounter. We can effec¬ 
tively vary the time since encounter by considering different 
angular locations in the Galaxy at the same radius as the 
Sun. The rotation period at the Sun is P ~ 0.23 Gyr, hence 
an uncertainty in the time since the El encounter of only 
0.1 Gyr would give rotation by nearly 180°. Estimates of 
the time since the last pericentre encounter range from 0.8 
to 1.1 Gyr. For example, Purcell et al. (2011) give a time 
and distance for last pericentre for the Sagittarius dwarf 
galaxy of t per i ~ 0.85 Gyr ago, R per i ~ 15 kpc for their 
light model and t per i ~ 1.1 Gyr ago , R per i ~ 15 kpc for 
their heavy model. Law X Majewski (2010) in their Figure 
7 give t per i ~ 0.8 Gyr ago, R per i ~ 15 kpc for the Sagittarius 
dwarf. 

In Figures 8, we emulate Figures 10, 12 and 13 by 
Williams et al. (2013) and Figure 2 by Carlin et al. (2013) by 
plotting gradients of (Vr) , (Vq) — V©,and (Vz) in different 
local neighbourhoods for encounter El in (R, Z) with 6 ^ 
R ^ 10, —2 ^ Z ^ 2 kpc and 0 in 45° increments with a 16° 
spread using bin sizes of (AR, AZ) = (0.25,0.25) kpc. We 
generate mean velocity matrices in the same method as pre¬ 
viously above. The resulting mean velocities are smoothed 
with a Gaussian filter with ay , &x = 1 to reduce the grain¬ 
iness resulting from the particle distributions. Like Figures 
4 and 5 we only show in Figure 8 the distribution at a par¬ 
ticular time since encounter, but the structure would be dif¬ 
ferent if the encounter occurred earlier or later. Figure 8 
shows that there is potentially quite a bit of structure in the 
bulk velocities above and below the Galactic plane, and that 
there could be large variations in the vertical gradients as a 
function of R , 0 as well as Z. If phase wrapping of epicyclic 
perturbations is responsible for the observed variations in 
bulk velocities, a consequence would be predicted variations 
in the bulk velocity as a function of distance from the Sun, 
such as the symmetry in (Vr) and (Vq) — V© above and 
below the midplane. 

Figure 8 illustrates that there are regions in our simu¬ 
lated galaxy that show negative vertical velocity gradients 
in (Vz). For example, the second panel corresponding to 
0 = 45°, displays negative Vz above the Galactic plane for 
R > 8 kpc and vice versa below the Galactic plane, giving 
Az with the expected and observed sign and about the right 
size. The velocities near the midplane in most regions are 
approximately zero, consistent with the presence of unper¬ 
turbed stars in the extracted neighbourhoods. The largest 
velocities are seen away from the Galactic plane, consistent 
with the presence of perturbed stars. However, as we did 
not include a thick disc component in our initial particle 
distributions, above and below the Galactic plane are found 
only perturbed particles and these dominate the mean ve¬ 
locities. Had we included a thick disc component we could 
have measured lower gradients as they could have diluted 
the contribution from the perturbed stars. We adopted a 
1/R number density to increase the number of particles at 
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Figure 7. Mock observations of the solar neighbourhood for breathing and bending mode parameters Az and Bz similar to Sun et al. 
(2015) Figures 13 and 14 for encounter El. Each subpanel is at a different central 0 and R with values for © and R shown on the left and 
bottom. The X axis of each subpanel has range —2 ^ Z ^ 2 kpc and the Y axis —30 ^ (Vz) ^ 30 km s _1 , the linear fits showing vertical 
gradients in Vz- This local neighbourhood shows smooth variations in mean velocity components as a function of Z. The gradients in 
the mean velocities vary rapidly with position in the simulated disc. 


large radius, and this too would bias our computed values 
of Az and Bz- 

Carlin et al. (2013) show negative (Vr) ~ — 5 km s -1 
near the midplane and positive {Vr) ~ 5 km s -1 above and 
below the midplane. Our integrations primarily show a zero 
value of (Vr) in the midplane because unperturbed disc par¬ 
ticles are in nearly circular motion. However, in many of the 
rows shown in Figure 8 the (Vr) distributions are symmet¬ 
rical at positive and negative Z , such as in the plot at 45°. 
We attribute this to the fact that v ~ 2k throughout the 
disc, so vertical epicyclic motions have sheared out more 
than have radial epicyclic motions, leaving the radial struc¬ 
ture more coherent. This symmetry is exhibited in numerous 
extracted regions in Figure 8, though the 45° panel shifts be¬ 
tween negative and positive Vr at about R = 8 kpc. Carlin 
et al. (2013) shows a negative value for Vr at \Z\ < 300 pc so 
the velocities could be associated with self-gravitating spiral 
density waves, and these are not present in our integrations. 
One last thing to point out is that the (Ve) — U© panels 
show regions that are both positive and negative above and 
below the plane (usually symmetrical above and below the 
plane). Asymmetric drift associated with a thick disc high 
velocity dispersion would lower the tangential velocity com¬ 
ponent. If perturbed populations of stars dominate at high 
and low galactic latitude, one would have to take care not 
to misinterpret the mean tangential velocity components. 

In summary, mock observations of our simulation at Rq 
and different azimuths show that perturbing the disc with 
a fairly massive Sagittarius dwarf yields measurable gradi¬ 
ents with similar sizes compared to recent observations. The 
vertical velocity component varies smoothly with Zi and has 
slope that varies quickly with position in the Galaxy (Figure 
7), mimicking measurements by Sun et al. (2015). Mean ve¬ 
locities as a function of (Z) (Figure 8) display gradients and 
symmetry above and below the midplane similar to those 
found by Carlin et al. (2013); Williams et al. (2013), how¬ 


ever we do not identify a particular region in our simulated 
galaxy that matches the gradients of all observed velocity 
components. 


4 SUMMARY AND CONCLUSION 

Non-interacting test particle integrations in a fixed poten¬ 
tial are less accurate than N-body simulations. However, 
they have the advantage that the role of phase-wrapping 
of epicyclic perturbations can be studied in isolation from 
phenomena that require self-gravity such as bending and 
breathing waves and modes. We start with a thin disc con¬ 
taining stars subjected to a velocity impulse computed from 
an instantaneous hyperbolic orbit approximation from a sin¬ 
gle close encounter with a dwarf galaxy. Particles are then 
integrated forward to the present time in a fixed Galactic 
potential and during this time the distribution of epicyclic 
angles is progressively sheared due to the dependence of ver¬ 
tical and radial epicyclic oscillation frequencies on mean or¬ 
bital radius. To illustrate this process we use two encoun¬ 
ters taken from estimated orbits for the Sagittarius dwarf 
galaxy nucleus. After the orbit integration the mean veloci¬ 
ties of the stars display vertical velocity gradients that have 
been previously interpreted in terms of breathing and bend¬ 
ing modes or waves. However, because our integrations lack 
self-gravity, these phenomena must instead be due to shear¬ 
ing or phase wrapping of epicyclic perturbations caused by 
the encounters. 

Because the encounter excites radial perturbations, the 
perturbed particle distribution not only shears with az¬ 
imuthal angle but wiggles, overlapping the distribution of 
unperturbed particles. We attribute the measured vertical 
gradients in bulk velocities to regions where perturbed and 
unperturbed particles overlap or are found at the same pro¬ 
jected A, Y in Galactocentric coordinates, or regions where 
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primarily perturbed stars are located, where underdensities 
would occur. A vertical gradient arises when a population of 
perturbed stars lies above or below the plane and its veloc¬ 
ity components are compared to those of unperturbed par¬ 
ticles that lie in the midplane. Vertical gradients could also 
arise when different populations of perturbed stars overlie 
the same X : Y location. This scenario predicts correlations 
between vertical gradients of different velocity components 
and density. The number of peaks in radial or vertical ve¬ 
locity means can be estimated from the time since the en¬ 
counter and the phases in epicyclic angles induced by the 
encounter. 

This work solely considers epicyclic phase wrapping on 
the Galaxy after an encounter from an external perturber. 
Effects from internal perturbations, such as spiral arms or a 
bar, or those due to self-gravity, such as breathing and bend¬ 
ing modes, are not present in our integrations. Because phase 
wrapping of epicyclic perturbations can produce vertical 
gradients in bulk velocities, breathing and bending modes 
are not required to account for these gradients. The wave¬ 
length of bending waves near the Sun has been estimated 
to be quite large (~ 10 kpc Weinberg 1991) implying that 
variations in bulk velocities over short distances must be due 
to coherent epicyclic motions, rather than self-gravitating 
waves or modes. That is not to say that self-gravity has no 
effect. Internal structures such as spiral density waves or 
the Galactic bar might also induce velocity gradients (Faure 
et al. 2014; Monari et al. 2015). An external perturber in 
the outer disc would produce gradients in bulk velocities 
that increase in magnitude with increasing distance from 
the centre of the galaxy, while an internal perturber (such as 
the bar) would exhibit gradients with decreasing magnitude 
with larger radius (e.g., Monari et al. 2015), so future ob¬ 
servations could differentiate between the two mechanisms. 
The radial distance in maxima of (Z), compared to bending 
wavelengths, may make it possible to distinguish bending 
waves from epicyclic phenomena. Unfortunately, studies of 
breathing modes have been restricted to the plane parallel 
setting (Widrow et al. 2014; Widrow V Bonner 2015) and 
so do not predict structure as a function of position in the 
Galaxy. 

We find that the simulated sizes of the Az and Bz coef¬ 
ficients, used to quantify gradients of bulk vertical motions 
by Widrow et al. (2014), are similar to those observed by Sun 
et al. (2015) for a moderate mass Sagittarius dwarf galaxy, 
~ 2.5 x 10 lo M Q , that passed its orbital pericentre about 
a Gyr ago. The passage of the Sagittarius galaxy through 
the Galactic plane approximately 2 Gyr ago is less effective 
than the pericentre at approximately 1 Gyr ago for three 
reasons. The extent of phase wrapping is more extreme and 
warping more tightly wound. Secondly, the E2 perturbation 
itself excites vertical epicyclic motions that vary in phase 
by 7r, whereas the pericentre encounter 1 Gyr ago, because 
it passes above the disc, excites vertical epicyclic motions 
in phase. Thirdly, even though we used a larger mass for 
the E2 encounter, the size of vertical velocity perturbations 
is smaller than that for the El encounter because of the 
orientation of the encounter. 

Regions at a solar neighbourhood radius extracted from 
our El integrations can show negative values of the Az 
coefficient, similar to that recently measured in the Solar 
neighbourhood (Carlin et al. 2013; Williams et al. 2013) but 


not at our expected location of the Sun in the simulation. 
However, a small error in the time estimated since the lat¬ 
est encounter could account for this discrepancy. We fail to 
find a local region in the El integration that matches the 
observed vertical gradients of all velocity components. Our 
integrations do illustrate that phase wrapping of epicyclic 
perturbations caused by the Sagittarius dwarf galaxy might 
account for much of structure seen in bulk velocity motions 
away from the Galactic plane, and if so, there should be large 
variations in the bulk motions with position in the Galaxy. 

We simulate two perturbations from the Sagittarius 
dwarf galaxy, each with different mass at different times 
and positions, from one orbit of the dwarf galaxy. If the 
dwarf is massive enough that dynamical friction is impor¬ 
tant, then the initial mass of the Sagittarius dwarf affects 
the orbit (Purcell et al. 2011), further contributing to the 
uncertainties in orbital parameters for our encounters. It is 
possible that the Galactic disc could have been more recently 
perturbed by an as yet unidentified dwarf, as proposed by 
Chakrabarti & Blitz (2009). 

In this paper we have used highly simplified integrations 
to isolate epicyclic phase wrapping from other phenomena. 
The impulse approximation for the encounters could over¬ 
estimate the energy transfer to the Galactic disc (D’Onghia 
et al. 2010), as disc response and time dependence of en¬ 
counters have not been taken into account. Our orbit inte¬ 
gration neglects spiral structure and associated radial mi¬ 
gration that would vary oscillation frequencies and so cause 
a loss of coherence in the epicyclic phases(e.g., Vera-Ciro 
& D’Onghi 2015). Self-gravity would also cause variations 
in the epicyclic amplitudes and angles and so might need 
to be taken into account to relate structure in the ve¬ 
locity field to a previous encounter. A future comparison 
between simulations with gravitationally interacting parti¬ 
cles and non-interacting particles could be used to identify 
and study possible bending waves or breathing modes that 
might be present in N-body simulations. Test particle sim¬ 
ulations could be used to improve understanding of the role 
of the time dependence of the encounters and relate the size 
scale and distribution of velocity gradients and bulk motions 
to the perturbations themselves. Thus, improved numerical 
models and observations might together match the observed 
gradients and simultaneously better constrain the mass and 
orbit of the Sagittarius dwarf galaxy. 
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APPENDIX A: INSTANTANEOUS 
HYPERBOLIC ORBIT APPROXIMATION FOR 
VELOCITY PERTURBATIONS 

We approximate the perturbation caused by the passage 
of the dwarf galaxy through the Galaxy midplane (E2 en¬ 
counter) and pericentre encounter (El encounter) with a 
velocity impulse. Each star’s velocity is changed instanta¬ 
neously. The velocity perturbation is estimated assuming 
that the encounter can be approximated with a hyperbolic 
orbit (e.g., see (Binney <V Tremaine 1987 section 7.1). We 
use an instantaneous hyperbolic orbit orbit approximation 
rather than the impulse approximation (Spitzer 1958) (the 
high velocity limit of the hyperbolic orbit approximation) 
because the E2 encounter velocity is oriented nearly per¬ 
pendicular to the disc. When the perturber velocity vec¬ 
tor is perpendicular to the disc, the impulse approximation 
gives no vertical velocity perturbation to stars in the disc, 
making it impossible to explore the role of vertical epicyle 
phase wrapping. However, the more accurate hyperbolic or¬ 
bit approximation gives a small velocity perturbation in the 
direction opposite to the velocity vector and so there is a 
vertical velocity perturbation. 

A low mass particle encountering mass M with a rela¬ 
tive velocity V and with impact parameter b on a hyperbolic 
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Figure 8. Mock observations of mean velocity components in local regions as a function of R and Z for simulation El. Each row shows 
(from left to right) (Vr), (Vz), and (Ve) — V c i rc (R). Mean velocities have been extracted at different azimuth angles. Panels from top to 
bottom have central angle © = 0,45,90 & 135°, respectively. The (Vz) panel at © = 45° displays a particularly large gradient. 
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Figure 8. Figure 8 continued, with (top to bottom) © = 180,225,270 & 315°. 
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orbit has, after the encounter, a change in velocity 


Av± 

Avw 


2 bV 6 


GM 

w(i + 


1 + 


b 2 V A 

(PM 2 

.2 t/ 4 x - 1 


6W 

G 2 M 2 


(Al) 

(A2) 


(Binney & Tremaine 1987 section 7.1) where Av± is in the 
direction toward position of closest approach and Avy is in 
the opposite direction of the relative velocity vector of M 
and the low mass particle. 

To compute the impact parameter and vertical and par¬ 
allel direction vectors we consider linear trajectories for both 
star and dwarf galaxy, x^(t) = x^o + Vd(£ — to) for the dwarf 
galaxy with position and velocity x^o and at time to when 
it passes through the Galactic plane. For a star at x s o and 
velocity v s a linear trajectory x s (t) = x s0 + v s (t — to). These 
linear trajectories are those in the vicinity of the encounter, 
ignoring the gravitational interaction during encounter and 
the background Galactic potential, and we only use them to 
compute the velocity perturbations. In the frame with rela¬ 
tive velocity V = — v s centered on the dwarf galaxy, the 

linear trajectory of the star is x S;Com = x s0 — x^o — V(£ — to)- 
By minimizing distance to the dwarf we find that the star 
(on the linear trajectory) is closest to the dwarf galaxy at 
time tmin = to + V • (x s0 — ^do)V ~ 2 . At this time the star 
is at 


X s,com{tmin ) — X s o X-dO V [V • (x s0 — Xdo)] U“ 2 . 

The vector between dwarf and star (and pointing toward 
the dwarf for the linear trajectory at closest approach) is 
the equivalent to this but with opposite direction 

b = x^o — x s0 + V [V • (x s0 - x d0 )] V~ 2 

The length of b is the impact parameter, b. The direction of 
b gives the perpendicular direction for the velocity change 
(equation Al). The parallel velocity perturbation (equation 
A2) is in the direction of the relative velocity vector —V with 
unit vector V = (v<j — v s )/|(vd — v s )|. Rewriting equations 
Al and A2, we can write the total velocity perturbation of 
the star as 

. 2bV 3 ~ /. b 2 V 4 V 1 

V "GM(6) b l 1 + Gnm*) 

- ( b 2 V 4 \ ' 

~ 2VV ( 1 + Gniw) 

with first term the perpendicular component and the second 
term the parallel component. Here b is the unit vector with 
the direction of b. We approximate the dwarf galaxy mass as 
that integrated out to the impact parameter, M(b) to take 
into account the spatial extent of the dwarf galaxy and to 
limit the size of the perturbations at small impact param¬ 
eters. A Hernquist mass model (mass interior to radius r) 


M{r) = M d y 2 (A4) 

(r + a H y 

is used to model mass distribution of the dwarf galaxy, with 
Hernquist scale length clh = 3 kpc. The scale length of the 
dwarf galaxy truncates the size of the largest velocity per¬ 
turbations. However, as large perturbations are only induced 


over a very small area, the velocities for most of the per¬ 
turbed disc are insensitive to the value of clh used. 

We restrict the effect of the hyperbolic approximation 
to stars within 18 kpc of the location of impact (the positions 
given in Table 1). Impact parameters larger than this (such 
as those for stars on the opposite side of the disc) would not 
be strongly perturbed by the encounter and the perturba¬ 
tions would not be well approximated by the instantaneous 
hyperbolic orbit approximation. 

In the limit of high velocity or b 2 V 4 / (GM (b)) 2 1 

equation A3 becomes 


A „ 2 GM(b) 

Av = v w (A5) 

resembling the impulse approximation for a point mass of 
mass M(b). For large impact parameter, b clh, the en¬ 
closed mass M(b) ~ Md and 


Av 


^ 2 GM d 


as would be expected for the impulse approximation for a 
point mass of mass Md. The impulse approximation (Spitzer 
1958; Cincotta et al. 1991) adopts a straight line for the the 
trajectory of the perturber (here the dwarf galaxy), with 
respect to a star in the disk. The velocity perturbation of 
the star caused by an extended potential perturber (rather 
than a point mass) can be computed as 


Av 



dtV$(t) 


(A6) 


(Cincotta et al. 1991) where <h(t) is the potential of the per¬ 
turber on a linear trajectory. This can be integrated analyt¬ 
ically for the Plummer model (Cincotta et al. 1991). Using 
the Hernquist gravitational potential, $(r) = —GMd/(r + 
clh), consistent with equation A4, equation A6 gives 

Av Kb ^^ g(a H /b) (A7) 

with function 

nUA= r _ dx _ 

9[U J 0 (\/PTT + u) 2 \/? r +T' 

In the limit of b clh, the limit lim u _>.o g(u) — 1 and we 
recover the impulse approximation for point mass M — Md, 
consistent with the total enclosed perturber mass at large 
radius. 

How accurate is equation A5 compared to equation A7? 
Computing the ratio of the two expressions at b — an we 
find that equation A5 underestimates |Av| by 30%. This 
underestimate is due to the neglect of the mass outside of 
radius r — bin equation A5. Drawing from this computation, 
we estimate that we underestimate the velocity perturbation 
at b — clh using the hyperbolic orbit approximation (equa¬ 
tion A3) by about the same fraction. For our two encounters 
most stars have impact parameters larger than clh so we do 
not expect this underestimate to significantly influence our 
results. 
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